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We model the quasicrystal-related structure CaCd6, a bcc packing of icosahedral clusters containing tetra- 
hedra which undergo orientational orderings at T<100 K. We use general schemes to evaluate an effective 
Hamltonian for inter-tetrahedron orientations, based on all-atom relaxations, either in terms of discrete cluster 
orientations, or of continuous rotation angles. The effective Hamiltonian is used in Monte Carlo simulations 
to find the (complex) ground state ordering pattern as a function of pressure. A preliminary investigation of 
thermal transitions found (in part of the pressure range) two different first-order transitions occurring below 100 
K. 

PACS numbers: 



I. INTRODUCTION 

Icosahedral i-CaCd (and isostmctural alloys such as i- 
CdYb) are the only known binary quasicrystals that have sta- 
ble long-range-order.——. Furthermore, their atomic structures 
represent a third family among quasicrystals, distinct from 
the previously known families of aluminum-transition metal 
and of Frank-Kasper packing (though including features of 
each). The structure contains large icosahedral "Tsai" clusters 
with tetrahedra of Cd atoms at their centers, which obviously 
breaks the cluster's symmetry, is built from from icosahe- 
dral "Tsai" clusters consisting of several concentric shells; the 
outer shells have icosahedral symmetry, but the innermost one 
is a Cd4 tetrahedron which can relatively easily rotate to dif- 
ferent orientations. This obviously breaks the cluster's sym- 
metry. 

This paper is concerned with the tetrahedra in CaCdg, a bcc 
packing of the same Tsai clusters; an equivalent phase is stable 
in many other systems that form binary quasicrystals (e.g. Cd- 
Y, Cd-Eu, etc). Periodic structures having a unit cell like a 
fragment of a quasicrystal phase are called "approximants"; 
CaCdg is the simplest approximant of the i-CaCd quasicrystal 
and others are known in nature. 

Diffraction on CaCdg-type approximants has revealed var- 
ious order/disorder transitions, which are ascribed to orienta- 
tional ordering of the clusters. Thus, as a function of pressure 
(up to 5 GPa), CdgYb has a complex phase diagram with six 
phasesi However, the exact orientations have not yet been 
determined experimentally or explained theoretically for any 
of these phases. In this paper, we compute a comprehensive 
set of interactions which, we hope, will predict the orienta- 
tional ordering patterns and transition temperatures, and may 
serve as a starting point to address the role of orientations in 
stabilizing the quasicrystal phase. 

The first modeling of tetrahedron energies in CaCdg — used 
ab-initio energies and started by assuming the experimentally 
refined sitesi^. Thus it was a sort of energy-guided fit, re- 
solving the correlations in partially-occupied sites found from 
diffraction. 

Later, Brommer et al performed a multiscale analysis to 



understand the interaction and show the transition behav- 
ior: first building inter-atomic potentials^ then modeling 
the nearest-neighbor interactions with 42 parameters (minus 6 
constraints)^, and finally performing some Monte Carlo sim- 
ulations with this effective cluster Hamiltonian-, Our paper 
should be considered a followup of this work. 

In this paper, we build a systematic method find the effec- 
tive interactions of tetrahedra (or similar inner clusters in other 
materials) from numerical relaxations. We give a natural way 
to eliminate arbitrary redundant freedoms in the interaction, 
so as to ensure the physical relevance of the fitted parame- 
ters in our model Hamiltonian, leading to a better view of the 
interactions (Section [IIP 3t . A singular- value decomposition 
is used to identify the dominant contributions in the potential 
(and, in principle, to reduce the number of terms needed to 
represent it). We also show (Sec. [IV} how to extend the same 
framework from a discrete set of orientations to the whole ori- 
entation space, and try to infer the functional form of the clus- 
ter orientation interaction Hamiltonian. Moreover, in Sec. IVl 
we use the effective Hamiltonian to find the lowest energy 
states in super cells not accessed by Brommer's original cal- 
culations^. We explore other low energy structures that are 
found in larger super cell Monte Carlo relaxations. 



II. FRAMEWORK AND METHODS 

In this section, we introduce the general concepts and meth- 
ods used in the rest of the paper, specifically the nature of 
the inter-cluster effective Hamiltonian (Sec. Ill Al l, the micro- 
scopic calculation of relaxed energies (Sec. Ill Bb . and our 
scheme for extracting and processing the effective Hamilto- 
nian (Sec. Ill Dt . Most of them are not specific to Tsai clusters 
with tetrahedra, but would also work for other kinds of reori- 
entable interior clusters inside icoahedral clusters, such as the 
pseudo-Mackay icosahedron in Allr or AlPdMn materials^. 

Before specializing to the clusters, we will review the con- 
stant structure that surrounds them. In CaCdg, it consists of a 
bcc packing of icosahedral Tsai clusters, with lattice constant 
15. 7 A. Each Tsai cluster consists of the following concentric 
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shells: 

(1) Z114 tetrahedron, radius 1.9 A; 

(2) Zn 2 o dodecahedral cage, radius ~ 4.2A; 

(3) Cai2 icosahedron, radius 5.56A; 

(4) Zri3o icosidodecahedron, radius ~ 6.4A. 

These clusters touch along the 3 -fold direction, while there are 
a few more Zn atoms between clusters around the 2-fold direc- 
tion. (Alternately, these Zn atoms may be reckoned to belong 
to large triacontrahedra of Zn on both vertices and midedges, 
which overlap along the 3-fold inter-cluster linkage.) 

A. Cluster effective Hamiltonian 

In this material, the low-energy degrees of freedom are 
the tetrahedron cluster orientations, represented by rotation 
matrices {Qi} relative to some reference orientation, where 
i = l...iV C eU as there is one cluster in each of the iV C eii prim- 
itive cells. The positions of all other atoms are taken to relax, 
so as to accommodate the tetrahedra^. they may, indeed, 
have large displacements, but these are dependent on {Oj}. 

We define an effective cluster Hamiltonian H({£2j}) as the 
minimum energy taken over all possible configurations con- 
strained to have that combination of orientations, allowing re- 
laxations of the surrounding atoms as well as distortions of the 
tetrahedra^i This effective Hamiltonian breaks up into one-, 
two-, and many-cluster terms: 

H(n) = n 1 + n 2 + n 3 + --- (2.1) 

We will assume that the many-body interactions are negligi- 
ble. Of course, the terms have the full symmetry of the crystal 
structure (minus the tetrahedra): for example, the one-body 
term is the same for all clusters and is invariant under the point 
group m3 (=T h ). 

There is a useful analogy between the cluster degree of free- 
dom and a (classical) spin, a unit vector that is specified by 
only two Euler angles, in contrast to the 3x3 rotation matrix £1 
which requires three Euler angles. The two-cluster interaction 
is analogous to dipolar or exchange spin interactions, while 
the single-body potential U(Q) is analogous to the single-spin 
spin anisotropies due to crystal fields. An even closer analogy 
is to the interacting, rotatable CN dipoles in KBri_ a; (CN) a ~2. 

1. Single-cluster terms and optimal orientations 
The single-body terms are 

i 

The single-cluster potential U(£f) includes a large contribu- 
tion with icosahedral symmetry, reflecting the strong steric in- 
teraction between the tetrahedron and the dodecahedral Cd2o 



cage around it. Because those cage atoms sit practically at the 
hard-core radius, it is not surprising if U(Cl) has some sharp 
and irregular-looking dependences on orientation. The single- 
body term is expected to have a somewhat smaller contribu- 
tion with cubic symmetry, indirectly due to the Tsai cluster's 
outer shell being distorted by its surroundings. 

These tetrahedra get from one orientation to another by a 
quasi-rigid rotation: rigid in the sense that the topological 
identity of the tetrahedron is maintained throughout, but in 
fact both the tetrahedron and its Cd2o cage undergo strong 
distortions. (Indeed, sterically the tetrahedron cannot never 
change orientations unless there are cooperative motions of 
the caging atoms.) A previous study— addressed the barriers 
and dynamics of a single Zn4 tetrahedron in the ZngSc ap- 
proximant (isostructural with CaCdg). The present paper is 
concerned only with static properties. 




FIG. 1: CcU tetrahedron in one of the twelve discrete optimal con- 
figurations. For visibility, the tetrahedron atoms are shown by large 
balls whereas those of the surrounding Cd2o cage are shown by small 
spheres. Bonds are drawn among the cage's atoms to highlight its do- 
decahedral shape. The orientation shown is +X r , relative to the axes 
indicated. 

Simulations 1 ^ found that tetrahedra in CaCdg tend to relax 
into one of twelve symmetry-related discrete orientations 
(for fi = 1 • • • 12), which must be minima of the single-cluster 
potential U(ft). The fact that similar orientations are seen, 
regardless of how neighboring clusters are oriented, indicates 
the single-body term is at least as strong as the two-cluster 
term. However, certain discrete orientations are unstable in 
the presence of certain backgrounds of uniform surrounding 
orientations (see Sec. Ill DTI below): so the single-body term 
is not much stronger than the cluster interaction. 

A tetrahedron in one of the ideal orientations is shown in 
Figure Q] One of its twofold (actually 4) symmetry axes is 
lined up with one of the cubic coordinate axes; in the figure, 
this is the axis coming out of the paper. Now, there is no 
4 symmetry in the icosahedron's point group (or more perti- 
nently, one in the 2 /m3 point group of the cluster center in 
cubic CaCdg). Hence, the two ends of that twofold axis are 
inequivalent. (Indeed, Figure[T]shows the front two atoms are 
lined up under two atoms of the cage, whereas the two back 
atoms are rotated 90°.) Hence, there are six possible direc- 
tions for that orientation axis which we label ±X, ±Y, ±Z. 

It is not quite stable for the two "front" tetrahedron atoms 
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to line up directly under the two nearby cage atoms. The 
cluster relaxes by rotating around the orientation axis by an 
angle of approximately ±15° - in either sense, thus sponta- 
neously breaking a symmetry and giving twelve symmetry- 
equivalent directions. We indicate the rightwards (clockwise) 
or leftwards (counter-clockwise) rotation, as viewed from the 
tetrahedron center, by a subscript r or I, so our complete la- 
bels a form are written +X r , etc. (Since orientations related 
by r H I are relatively close, we anticipate they may have 
similar interactions.) 

In earlier experiments, Gomez and Lidin studied the x-ray 
diffraction of MCd 6 approximants, where M= Ca, Y, or rare- 
earth. They mapped out the continuous electron density inside 
Tsai clusters, which they were able to interpret in terms of a 
host of split positions representing tetrahedron orientational 
disorder, with preferred orientations of a single kind related by 
symmetry^. The apparent symmetry (see Ref.[jj| Fig. 3) de- 
pended on the kind of large ion M, presumably reflecting the 
relative importance of the icosahedral and cubic components 
in the single-body term U(SY) of the orientational Hamilto- 
nian: icosahedral for M=Tm and Lu, cubic for M=Tb, and 
intermediate for M=Ho or Er. Their result for the case CaCdg 
agrees with the simulations of Ref. [l6| as confirmed by our 
own: the tetrahedra sit in an asymmetric orientation. 



2. Cluster pair terms 
The two-body term is written 

« 2 =x;vy(a,a). (2.3) 

id 

The function is translationally invariant, depending on 
sites (i, j) only through the vector connecting them. It is ex- 
pected to decay with separation. In this paper, the only separa- 
tions included in the fit are the two kinds of nearest neighbors 
in the bcc lattice of cluster centers: the "6" linkage (separation 
vector equivalent to [0, 0, 1]) and the "c" type linkage (separa- 
tion vector equivalent to [1/2, 1/2, 1/2]). 

We now comment on the possible atomic-scale origins of 
the cluster effective interaction; however the results of this 
paper do not depend on understanding that, nor will they re- 
solve it. A priori, one expects the cluster pair interaction has 
two kinds of contributions: mediated elastically, via displace- 
ments of intervening atoms, or mediated by the electron sea. 
The latter is expressed, within our framework, by the EAM 
or pair potentials (see Sec. Ill B II below) and more specifi- 
cally by the long range Friedel oscillations characteristic of 
pair potentials in a metal. In contrast to the Al-TM and F- 
K classes of quasicrystal, Friedel oscillations do not appear 
to be crucial for the "Tsai cluster" class of quasicrystals 15 , 
which suggests the elastically mediated interaction should be 
dominant. Furthermore, if direct electronic interactions were 
dominant, one would expect the interaction of two clusters 
separated by a [1,0,0] type bond to be invariant under a si- 
multaneous rotation by 90° around the bond direction, which 
is not the case (see Sec. Hill ). However, ab-initio calculations 



by Ishii and collaborators ound a substantial cluster interac- 
tion when atom positions are not relaxed, demonstrating that 
the direct electronic interactions are significant. Finally, we 
remark that in the isostructural compound ScZng, the cluster- 
cluster interaction is much smaller—. 

Since the two-cluster term is mediated by a comparatively 
small distortion of the outer shells of the Tsai cluster, and/or is 
a sum of several potential terms, we anticipate that it is more 
smoothly behaved, and that the contributions from different 
neighbors will be additive. 

B. Interatomic potentials 

In order to do our fitting, we must build a database of re- 
laxed energies coming from a lower, more exact level of de- 
scription. We defined the cluster Hamiltonian as the relaxed 
minimum energy of the system, having fixed the orientation 
of every cluster. This opens up three questions: 

(1) How do we define or compute the energy of an arbi- 
trary atomic configuration? (Sometimes these energies 
can be computed directly from ab-initio relaxations, but 
here we needed to use "classical" potentials. 35 ) 

(2) considering that the tetrahedra are typically distorted, 
what is our precise definition of cluster orientations? 
(This is required in order to define the family of con- 
figurations we are minimizing over.) 

(3) How do we implement this constrained minimization 
reliably? 

This section explains our answers to each question including 
the important technicalities. 

1. Potentials 

Classical potentials are essential in various situations when 
molecular dynamics or relaxation is required in supercells 
containing many clusters. For the present problem, we used 
the minimum possible supercell which is 3 x 3 x 3 or about 
4000 atoms, which is too large for doing repetitive ab-initio 
relaxations. 

Most of our calculations are based on the embedded-atom 
method (EAM) potentials fitted by Brommer and Gahler— . 
Their method of fitting is laborious and cannot be quickly 
repeated for a new material. As an alternative, we also 
tried the empirical oscillating pair (EOP) potentials 17 , which 
can be rapidly computed for any composition, but are valid 
only while the conduction electron concentration is held con- 
stant. Comparing results from the two potentials, as done in 
Sec lIIICI below, may give a measure of the uncertainty of our 
conclusions, and/or a measure of the reliability of plain pair 
potentials in this system where their validity is less assured. 

The (EOPP) pair potentials use an six-parameter analytic 
form which incorporates Friedel oscillations-!" 7 .. For this pa- 
per, they were fitted against a database of ab-initio results 
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FIG. 2: Fitted "empirical oscillating pair potentials" (EOPP) for the 
Ca-Cd system. 



with a total of 28 energy datapoints taken from relaxed 
T = OK structures, plus a single snapshot from a high- 
temperature molecular dynamics simulation at 1000K (of the 
cubic CaCdg structure) that gave over 7000 force datapoints. 
Our database of relaxed samples included all known Ca-Cd 
binary compounds (CaCd2 in both CeCu2 and MgZn2 struc- 
tures; Ca3Cd2, Ca2Cd7, and the B2 structure of CaCd; and 
versions of the CaCd^ approximant with six different ways 
of placing two Cd4 tetrahedra in the cubic 1/1 cell. Further- 
more, we added structures that we took from similar systems, 
such as Ca3Zn, Ca2Cu, CaZn3, CaCdn and finally the Frank- 
Kasper "Bergman" phase strucure of AlMgZn. In the final 
iteration, the fit was biased so as to give forces as accurate as 
possible. The results are shown in Figure |2] 

It should be noted that the Tsai class of quasicrystals (and 
related alloys) is based on either Cd or Zn, both of which are 
known in their elemental forms to have an anomalous c/a ra- 
tio in the hep lattice. Indeed, this might be related to the ex- 
tremely non-rigid behavior of the Zn2o or Cd2o dodecahedral 
cages in the Tsai clusters, which is essential in allowing the 
inner tetrahedra to rotate at all. (See frames from the finite- 
temperature MD simulation of ScZn 6 , Figure 1 of RefJ^). We 
also know that pure Zn is one of the few cases in which the 
EOP potentials more or less fail 17 . Consequently, it is some- 
what surprising that we find the EOP potentials succeed in 
CaCd6, in that the cluster Hamiltonian fit is similar to the re- 
sult from EAM potentials (see Sec. IIIICI ). 



1. Defining cluster orientations 

It is relatively easy to define the orientation O: the four 
inner atoms always do form some kind of tetrahedron, since 
it is sterically impossible for one of the atoms to pass through 
the plane formed by the other three. 

Let ri be the position of tetrahedron atom I (for I = 1 , . . . , 4), 
and define the center as r c = (fx + ?2 + f*3 + f^/A; note that 
f c can deviate from the center of the surrounding cage. Also, 
define a regular reference tetrahedron by four unit vectors {ii } 
in tetrahedral directions. (In the relaxation code, each tetrahe- 
dron's initial prescribed orientation is used as the reference.) 

Define a matrix M_ with components 

3 4 

M a0 = - Y^{?l ~ r c ) a {ti)p. (2.4) 
i=i 

Now write the polar (=singular value) decomposition 

M = Q.L Md Q-r (2-5) 

where Vl L and fl R are 3 x 3 rotation matrices, and M_ D is 
(positive) diagonal. It is easy to check that if {fi} is a regu- 
lar tetrahedron, then M D is a multiple of the identity and the 
actual tetrahedron is rotated, relative to the reference tetrahe- 
dron, by 

O = Q L Q R (2.6) 

For a general tetrahedron, we take Eq. (12.6) as our definition of 
its orientation. It can be shown that ( 12.6b optimizes a measure 
of agreement, ^? =1 (rj — f c ) ■ f2t/. If we had to do this for 
some other cluster with n a atoms, we simply need a different 
set of ideal vectors and the above recipe still works, with the 
replacement 3/4 — > 3/?i a in 12.41 . 

There is an alternative way to think of orientation extrac- 
tion, which is specific to the (present) case that the cluster has 
exactly four atoms, not in the same plane. We can uniquely 
and exactly represent the actual coordinates as 

n = Mti + r c (2.7) 

since this is a set of 4(3) linear equations in 3 2 + 3 unknowns 
(the components of M and r c . In materials science, such a 
matrix defining an affine transformation of the atoms is called 
the deformation matrix. Indeed, for four atoms it gives the 
same result for M_ as ( 12.4b . As above, a polar decomposition 
( 12. 5i is performed to define ft. 



C. Implementing cluster orientations and constrained 
relaxation 

We need to establish an explicit practical mapping between 
atom positions and cluster orientations, the degrees of free- 
dom at the two levels of description we want to relate. First we 
lay out the atoms-to-orientation mapping, and then its inverse, 
which actually means defining constraints for relaxation. 



2. Constrained relaxation 

In each constrained relaxation iteration, we assume a start- 
ing configuration in which each cluster already has its pre- 
scribed orientation, according to the definition based on polar 
decomposition. Conjugate gradient directions are constructed 
in the standard (unconstrained) fashion, but are then projected 
orthogonally into the allowed subspace, and one-dimensional 
minimizations are carried out along this projected direction. 
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We relax atom configurations, subject to the EAM or 
pair potentials, using a nonlinear conjugate gradient algo- 
rithm with Newton-Raphson and Fletcher-Reeves^. Each 
successive conjugate-gradient iteration consists of a one- 
dimensional Newton-Raphson minimization (with up to 10 
iterations) along the next conjugate-gradient direction. The 
stopping criterion is that the energy change per step is 
AE'stop < 10 _4 cV. Typically, a few hundred conjugate- 
gradient iterations were needed. 

The question is how to constrain all tetrahedron orientations 
while relaxing all atom coordinates. This amounts to 
3iV ce u nonlinear constraints, defined implicitly by (12.5) and 
(12.6b . The basic approach is to linearize these constraints, 
defining a linear subspace within the manifold of all atomic 
coordinates. 

Since the actual constraint is nonlinear, after some itera- 
tions the configuration would not exactly satisfy it. There- 
fore, every ~ 10 iterations we perform a nonlinear projection 
to reassert the orientation constraint. Namely, the actual re- 
laxed positions {r[} of the atoms in a tetrahedron are related 
to the ideal rotated positions {Oi/} by a deformation matrix 
M ', then polar decomposed as M_' = il L AI D il B as in ( 12. 5b . 
We then replace AP_ -> M_ D 

Specifically, let M be a possible additional infinitesimal 
deformation of the four atom positions {?^} after a step, rela- 
tive to the previous positions {r}}. 

f - f c = M'(f - r c ), (2.8) 

with A£ — I + m' with m' small. Then the condition that M ' 
contains no rotation is that mf is symmetric, which implicitly 
defines a set of linear conditions on f — r c . 



3. Comparison to unconstrained relaxation method 

We want to constrast our constrained relaxation app roach 
with the simpler alternative approach used in Ref. [l6f plain 
unconstrained relaxation starting with prescribed initial orien- 
tations (preferably, optimal ones). One difference is that our 
constrained relaxation allows construction of a full Hamilto- 
nian as a continuous function of orientations, as carried out 
below in Sectionsec:results-continuous, which allows simu- 
lating the orientation fluctuations found at T > or mapping 
out the energy barriers for a cluster to reorient. 

Using unconstrained relaxation, one accesses only the dis- 
crete set of orientations that are local minima. Usually, one 
implicitly depends on the assumption of a one-to-one corre- 
spondence between the combinations of initial orientations 
and final ones. In reality, it may happen - and did so in 
our study, in a few instances - that certain combinations of 
the discrete orientations in neighboring clusters are not lo- 
cally stable: they relax into some other combination of ori- 
entations. A final, technical drawback to unconstrained relax- 
ation is that even if that one-to-one correspondence exists, the 
actual ground state orientations deviate from the single-body 
optimum ones. 



D. Extracting the cluster Hamiltonian 

It is assumed we start by choosing a discrete list of repre- 
sentative possible orientations {u a , a = l...m}. In this work, 
these are either the twelve optimal orientations of the single- 
body terms (to obtain the discrete Hamiltonian in Sec. ITTTb or 
else a finer-spaced grid that is meant to sample a continuous 
range of orientations (in Sec. [IV}. Except where noted, our set 
of orientations always has the full point symmetry of 

cluster site. 

A dataset is then constructed of the relaxed energies E a p 
for orientation u a in cluster site i and in cluster site j, for 
all combinations of the two, while the other ("background") 
clusters are held fixed. In CaCd 6 , we consider (mainly) two 
kinds of pairs, which are the nearest-neighbor sites of the bcc 
lattice: those separated by (100) vectors and by (1/21/21/2) 
vectors. (We sometimes call these, respectively, "b" and "c" 
linkages, based on their role in the canonical-cell tiling^.) 

Our aim in fitting is to convert the array E a p to an energy 
function in the form of the cluster Hamiltonian (12. U . with 
well-defined and compact formulas for the single-body term 
U(P) and pair term V l3 (fi, Q'). 



1. Uniform background approach to isolate pairs 

In order to properly fit the pair interaction of two clusters, 
they must be put in a sufficiently large supercell that they have 
only one significant interaction with each other (no interac- 
tions with an image of the neighbor through periodic bound- 
ary conditions in a different direction). When the cluster pair 
is related by the b type vector [0,0, 1 ], that demands a supercell 
of at least 3x3x4 basic cells (i.e. 72 clusters). 

It is not feasible to exhaustively enumerate all possible 
combinations of cluster orientations in the supercell and re- 
lax every configuration; our fits must be based on a subset of 
orientations. We will describe two different ways to choose 
this subset: exhaustive enumeration of a pair in a fixed back- 
ground, which we used in this project, or random configura- 
tions of the whole system, used by Brommer et al^. 

Each supercell configuration is specified by three orienta- 
tions: fijjfij-, the orientations of the two clusters whose in- 
teraction we want, and the background orientation 0^, which 
is taken by all other clusters. While keeping £2 6 fixed, we 
enumerate all m 2 possible combinations of (0^ , ) and (as 
described shortly below) extract a fit V a p. Thus, we get m in- 
dependent fits, one for each background. (Due to symmetries 
not all of them are independent.) This provides some useful 
checks. 

When our sampled orientations {£2 Q } are just the twelve 
optimal ones, which are symmetry equivalent, U(£l a ) has 
the same value for every one so the correct result should be 
U a = 0. However, the presence of a particular background 
breaks the symmetry; in our fit, the pair interactions between 
cluster i or j and "background" clusters all get absorbed into 
U a term, so it will be non-constant. But if we average U a over 
all possible backgrounds, the symmetry should be restored. 
After we complete the fit for V a p, we can predict the apparent 
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single-body term due to the background, and check if this is 
consistent with the apparent single-body terms that were ac- 
tually fitted. (Any disagreement suggests the importance of 
further-neighbor interactions with the background.) 

In any case, provided we only have pair cluster interactions, 
we should still obtain an identical fit for V a p regardless of the 
background. But actually, multi-body interactions including 
both of the selected clusters, and one or more of the back- 
ground clusters, many contribute to the fitted V a p- Thus, any 
dependence of V a /3 on the background signals the presence of 
multi-body interactions. 

The uniform-background procedure relies on being able 
to generate any specified combination of orientations. This 
might fail, if we use unconstrained relaxation, because cer- 
tain combinations of orientations might be unstable and turn 
in a different direction. In particular, we found that the back- 
ground orientation £2 fcg is +X r , a varied cluster with orien- 
tation — X r very often is unstable, reorienting to —Xi. This 
problem is fixed by using a relaxation constrained to a particu- 
lar orientation (see Sec. Ill C 2l above). Alternatively, it may be 
evaded by using the random-whole-system approach instead. 

The random-whole-system approach was used by Brommer 
et al£ for CaCd 6 . (It was also used— in Ir 2 3Al 4 , a quasicrys- 
tal approximant from the Al-transition metal class, which has 
a different kind of orientable inner shell consisting of Aliolr 
in the pseudo-Mackay icosahedron cluster.) Produce 10 3 -10 4 
realizations of a supercell with a random combination of ori- 
entations and relax each one; if orientations are unstable, it 
is necessary to analyze the final state to ascertain the actual 
index a for each cluster. In each configuration, we find the 
count N a g of cluster linkages of type I having the clusters in 
orientations £2 Q and £2g, respectively, Then we use a linear 
least-squares fit to express the total energy as 

l,a,P 

Of course, symmetries are normally used so as to reduce the 
number of parameters. 

We mention briefly a third way to construct a database, re- 
lated in spirit to the uniform-background. (It is analogous to 
determining spin-spin interactions in a magnetic material by 
diluting it with nonmagnetic ions that are chemically equiva- 
lent to the magnetic ones, so that only isolated pairs of the lat- 
ter occur.) We choose two clusters to vary through all combi- 
nations of interactions, but take a background having no orien- 
tational interactions, by replacing every other tetrahedron by a 
single large atom that takes up the same space; as it happens, 
the Ca atom is close to the right size. There is no averaging 
over backgrounds since there is only one kind. 

The single-atom replacement trick was successfully applied 
to the isostructural compound ScZng, but only to study the 
effects of the single-cluster term on the dynamics (the inter- 
cluster interactions are, in any case, weak in that material). We 
experimented with this method for pair interactions in CaCd 6 , 
but these preliminary investigations suggested that the system- 
atic errors are too large for us to trust the results quantitatively. 



2. Symmetry properties of the discrete interaction matrix 

Consider the permutation symmetries of the interaction ma- 
trix V_ that follow from the space-group symmetry operations 
of the CaCdg framework, which affect V_ in two ways. First, 
the rotation part of these operations permutes orientations 
within the set {£2 a } (except in the trivial case of pure transla- 
tions). Second, there are three interesting cases for the action 
on the cluster centers: both may map to themselves, they may 
be swapped, or one maps to itself and the other maps to a dif- 
ferent cluster. For each of these cases, it will be convenient 
to write E a p and V a p using matrix notation as arrays E_ and 
V_, which in general are not symmetric. The action of a point 
group operation g on cluster orientations is written with an 
m x m permutation matrix, P_ g . (Here we mean the matrix 
in which each row and column contains exactly one element 
1, and the others are zero, so multiplication by this matrix in- 
duces a permutation of the indices.) 

First, for a lattice symmetry operation that keeps both clus- 
ters fixed, invariance of the energy under this symmetry is rep- 
resented by 

P g V(E g ) T = V. (2.10) 

This requires that the bond vector is invariant when g acts at 
the midpoint. 

Second, consider a rotational symmetry h that exchanges 
the two clusters. (In practice we need only consider inversion; 
all others are products of inversion and a symmetry of the first 
kind.) Now we get 

P h V{P h ) T = V T . (2.11) 

Thirdly, consider a rotation g that, acting on site i and maps 
i—ti but sends j —> k. We get 

P g V(P g ) T ^V (2.12) 

The use of ( 12.121 ) is that, having computed the interaction ma- 
trix for a pair separated by (say) [0, 0, 1], we find the interac- 
tion matrix for symmetry related bond vectors such as [100]. 



3. Redundant parameters: pair interaction resolution 

Whenever one fits an effective Hamiltonian to configura- 
tions of discrete variables, one finds linear dependencies be- 
tween the counts of certain local patterns and and certain other 
ones. (Here "local pattern" means a particular combination of 
two or more discrete variables in nearby sites.) Thus, if every 
local pattern is allotted an independent coupling coefficient, 
the values of these coefficients will be ill-determined: a whole 
family of parameter sets gives identical energies for all possi- 
ble configurations. Such dependencies arise when the discrete 
variables are Ising spins^ or tiles ("decorated" by atoms) in 
quasicrystal-related phases^; they arose in the previous fit of 
a cluster Hamiltonian 1 ^ and were handled by arbitarily setting 
certain coefficients to zero. 
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In fact, interacting clusters are the simplest of the cases 
where these dependencies are observed. Our task is to write 



E a p — E + U a + U' 3 + V a p. 



(2.13) 



The difficulty is that the solutions of (12.13b are ill-defined, 
since it is invariant under 



V a @ — > V a jj — A a — Ap] 

U a — > U a + A a : 



U' 



K + A' 



(2.14a) 
(2.14b) 
(2.14c) 



where {A a } is an arbitrary set of energy shifts. Within an ex- 
tended crystal, with bonds oriented in several directions, there 
will be different functions A a and A' a associated with each 
direction of bond; the one-body term gets shifted by the sum 
of all these. Thus, even if we constrained all U a = 0, that 
would not resolve the indeterminacy in { V a p} There is a sim- 
ilar indeterminacy in the single-cluster term: 



U a -B; 



E Q E Q + B + B'. 
We do resolve it by imposing the simple rule 



\v afi = Y i v a 



0; 



E^ - o- 

p 



(2.15a) 
(2.15b) 
(2.15c) 



(2.16a) 
(2.16b) 



Given an initial set E a p, it is easy to implement Eq. ( 12. 16ab 
using Eqs. ( 12.14b with 



= —(y^Eap); 
m — ' 



-4', 



(2.17a) 
(2.17b) 



After that we enforce Eq. (12.16b) by applying (12.15) with B 
QfL U a )/m and similarly B' . 



E. Singular Value Decomposition methods 

The matrix of interactions V a p for the discrete orientations 
has m 2 entries and this implies a large number of fit parame- 
ters, even after symmetries are taken into account and redun- 
dancies eliminated. For example, Ref. [l6| - taking the dis- 
crete approach - fitted a separate interaction parameter for ev- 
ery inequivalent combination of orientations of the nearest or 
second-nearest cluster neighbors, which amounted to 42 inter- 
actions (slightly fewer after removing redundancies). 

However, we anticipate that some combinations of these 
parameters are unimportant. Imagine, e.g., if a cluster de- 
veloped a charge q a dependent on its orientation, then the 
cluster orientation would be a Coulomb interaction propor- 
tional to q a qp/\Rij\ and thus of completely separable form. 



To parametrize such an interaction for m orientations, we do 
not need 0(m 2 ) numbers, but only one overall strength and a 
list of 0(m) charge strengths q a . An electric dipole interac- 
tion is similar, but since the dipole moment has three compo- 
nents, the interaction would take the form 

Our recipe to reduce parameters is to posit that the interac- 
tion matrices can be put in the form 



V, 



a/3 



£ G ^a9p 



(2.18) 



The implicit idea here is that the clusters are weakly perturb- 
ing some "field" which fills the space between them. Think of 
g£ as the "charge" presented by the first cluster when it is in 
orientation a, similarly g'Jf as the "charge" presented by the 
second cluster, while cr M is the effective interaction of these 
charges (expected to decay with separation). More exactly, fi 
indexes different kinds of "charge", or different components 
of the same charge. For example, if the tetrahedron atoms car- 
ried a literal charge, then /i would index the cluster's moments 
(monopole, dipole, etc) and their tensor components (three for 
the dipole case). In CaCdg, we expect the interaction is medi- 
ated by elastic strains of the intervening atoms, or possibly by 
direct couplings (which are represented by the oscillating tail 
in the pair potentials shown in Figure [2] or analogous tails in 
the EAM potentials^). 

For definiteness, let's augment ( 12.181 1, by an orthonormal 



condition J2 a 9 nag. 



up 



J flU- 



Then, these two equations are 
equivalent to the singular value decomposition, or in matrix 
notation 



V 



T I 

9 <L9 



(2.19) 



where a_ is the diagonal matrix g and g' are orthogonal matri- 
ces of singular values (which we take to be nonnegative and 
decreasing). In this form, we have m' = m& But our expec- 
tation (to be checked!) is that the singular values rapidly 
get smaller, so we can truncate (12.18b after (say) three terms 
and reduce the 42 parameters to a handful. 

In practice this decomposition is much more powerful when 
symmetries come into play. Following the analogies to (elec- 
tric or elastic) dipoles and quadrupoles, we imagine that the 
dependence of the "charge" on orientation il could be written 
as a smooth function g(d) with the angular dependence of an 
orientational harmonic. This idea is tested in Section HVI 



III. RESULTS: DISCRETE ORIENTATIONS 

In this section, we use the framework of Section Q]] to fit a 
cluster potential for the twelve discrete, symmetry -related op- 
timal orientations, similar to the form used in Ref. [HI Our 
microscopic Hamiltonian is defined by EAM potentials of 
Ref. [3 For each relaxation, we use an initial structure in 
which the "regular" atoms (all atoms excluding the tetrahe- 
dra) are in their averaged positions with cubic crystal symme- 
try (space group Im3), as refined in Ref. [US (this structure has 



8 





+X r 




+x t 


-Xi 


+x T 


-2.8 


0.5 


4.2 


-0.4 


+ Y r 


-8.6 


5.7 


-1.9 


3.9 


+ Z r 


-5.2 


2.0 


-1.1 


3.9 


— ■ JC r 


-6.8 


3.9 


-7.3 


4.2 


-Yr 


6.0 


-6.3 


3.2 


-1.1 


— Zf 


0.3 


2.6 


3.2 


-1.9 


+x t 


-5.3 


3.5 


3.9 


0.5 




-2.3 


11.7 


2.6 


2.0 


+z t 


-17.6 


11.7 


-6.3 


5.7 


-Xi 


9.6 


-5.3 


-6.8 


-2.8 


-Y, 


21.3 


-17.6 


0.3 


-5.2 


-Zi 


21.3 


-12.3 


6.0 


-8.6 



TABLE I: Interaction matrix for two clusters separated by 
[1/2,1/2,1 /2] (=c-linkage), using EAM potentials (in meV). Due to 
the threefold symmetry of this linkange, all interactions are invariant 
under cyclic permutations (xyz) in the labels when applied to both 
clusters; thus columns are given here only for the "X" orientations. 



the Pearson symbol cI208). Each tetrahedron is initially reg- 
ular (having radius 1.775 A and centered on the Tsai cluster 
center) and in the prescribed orientation f^. 



A. Cluster pair interaction 

We computed interactions for two kinds of cluster pairs, 
separated by the vector [1/2,1/2,1/2] (=c linkage) or [0.0,1] 
{=b linkage), placed in a 3x3x3 supercell. We used the uni- 
form background method and averaged over the twelve pos- 
sible backgrounds; as requiared by symmetry, the averaged 
single-body term was zero (i.e. all orientations were equiva- 
lent). The resulting pair interactions are given in Tables [F]and 

El 

To make sense of the computed interactions, it is essential 
to consider the local symmetry of the cluster pair, which is 
reflected in the symmetry of V a p as laid out in Sec. Ill D 21 

For the c pair, the point group (centered on the linkage's 
midpoint) 3 (-D3); the subgroup which maps each cluster to it- 
self is 3 (C3). This threefold rotation, in terms of our orienta- 
tion labels, just cyclically permutes (XYZ) without changing 
the ± or r/l part of the labels. As for the b pair of clusters, its 
point group is 2 /mm (£*2h) and the subgroup that maps each 
cluster to itself is 2m {p2v) which has order 4. Its generators 
are the mirror operations m x and m y , which respectively in- 
vert the x or the y coordinates. (Remember we singled out the 
z direction by aligning the pair with it.) The action of m x , for 
example, is to switch +X <R- —X, not switch the sign of the 
Y or Z orientations, and in all cases to switch r <-> I. Finally, 
in all cases it is convenient to let inversion be the fundamental 
operation that swaps the two clusters; its action on our labels is 

to always switch + < and always switch r -<-> I. (Note that 

any proper rotation always preserves the r/l indices whereas 
any improper rotation always switches them.) We take ad- 
vantage of these symmetries to reduce the number of columns 
displayed in the interaction matrices. 





+X r 


-Xi 


+ Y r 


-Y 


+Z r 


-Zi 


+X r 


17.4 


-17.0 


6.7 


6.5 


-1.1 


0.8 


+Xi 


18.7 


-16.5 


6.5 


6.7 


-0.6 


0.4 


-Xi 


-17.0 


17.4 


-7.0 


-8.5 


-0.6 


0.4 


— X r 


-16.5 


18.7 


-8.5 


-7.0 


-1.1 


0.8 


+ Y r 


-8.5 


6.5 


-2.5 


-2.3 


0.6 


-0.7 


+Y, 


6.5 


-8.5 


4.4 


2.8 


1.4 


-1.4 


-Y r 


7.0 


6.7 


-2.3 


-2.5 


1.4 


-1.4 


~Y r 


6.7 


-7.0 


2.8 


4.4 


0.6 


-0.7 


+Z r 


0.4 


0.8 


-1.4 


-0.7 


-3.8 


4.1 


+Z l 


0. 


0.4 


-0.7 


-1.4 


-2.6 


4.0 


-Zi 


-0.6 


-1.1 


1.4 


0.6 


2.6 


-3.8 


-Z r 


-1.1 


-0.6 


0.6 


1.4 


3.0 


-2.6 



TABLE II: Interaction matrix V a p for two clusters related by [0, 0, 1] 
vector (=6-linkage), in meV. In all interaction matrices the cluster 1 
orientation is given in column headings, and cluster 2 orientation by 
row headings. The interaction is invariant if one switches r-f+1 in both 
labels, so the columns with and "— r" labels are omitted. 



1. Justification of functional form 

Can we make sense of the patterns of interactions? First 
of all, we see that just a few combinations have much larger 
interactions than any others. In the c-linkage case, for ex- 
ample, the pair (— Zj, +X r ) and permutations has interaction 
21.3 meV; due to inversion symmetry relating the two clus- 
ters [we are using (12.1 1H . this has the same interaction as the 
pair (—Xi, +Z r ). In turn that is equivalent by cyclic permu- 
tation to (— Yi, +X r ) in our table, which is indeed 21.3 meV. 
The interaction (— Yi, —X r ) is almost as big, —17.6; this is 
equivalent to The interaction (+Xi,+Y r ), i.e. (+Zi, +X r ) 
in the table, by inversion. Overall, it can be seen that the r 
orientations have stronger interactions than the I orientations. 

Similarly, in the b interactions (Tablelllli. we that X orienta- 
tions (of all flavors) have by far the largest interactions, while 
Z orientations have the least. 

If the clusters were far apart and weakly perturbing their 
surroundings, we would expect the inversion of a cluster 
would have exactly the opposite coupling, in the case of a 
tetrahedron, since it moves the atoms to the places comple- 
mentary to them. To say this another way, consider the union 
of a tetrahedron and its inverse, i.e. a cube: whatever is the 
lowest nonzero "multipole moment" of the tetrahedron, we 
would expect that one to be zero in the case of the cube. Thus, 
we expect that approximately V a p = —V a >p when a and /3 
are related by inversion. This expectation is borne out in the 
6-bond interactions, as can be seen by comparing adjacent en- 
tries in columns 1 and 2, 3 and 4, or 5 and 6. However, in 
the case of the c-linkage it does not work well: the largest 
interactions ((—17.6 and +21.3) are related (see the leftmost 
column) by +Zi <-> — Z\. 

A somewhat surprising fact is that the size of c-bond in- 
teractions is the same as that of the fe-bond interactions; this 
is confirmed (see below) by the fact the respective dominant 
singular values from each interaction are practically the same. 
The best absolute measures, the matrix norms, are respec- 
tively 94.7 and 83.1, for a ratio ~ 1.1. One would have 
expected the c-bond interaction to be larger - by a factor of 
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a (meV) 


irrep 


+X r 


-Xi 


+x l 


— X r 


81.51 


A 


-0.1380 


-0.0054 


-0.3181 


0.4616 


33.79 


E 


0.2339 


0.4564 


0.4140 


0.4819 






-11.7° 


164.9° 


-66.9° 


138.0° 


5.443 


A 


0.4204 


0.0739 


-0.3675 


-0.1268 


2.278 


E 


0.7060 


0.1414 


0.2760 


0.2684 






177.1° 


-68.0° 


-109.8° 


161.6° 


1.385 


E 


0.2987 


0.2812 


0.5896 


0.3883 






142.7° 


176.5° 


-68.5° 


-79.5° 


1.330 


A 


-0.2329 


0.4944 


-0.1170 


-0.1445 


3.288 x 10" 2 


E 


0.1557 


0.5994 


0.2673 


0.4600 






33.4° 


65.8° 


55.8° 


-56.0° 



TABLE III: Singular values and right singular vectors for c-linkage. 
This cluster pair has a symmetry under cyclic permutations (XYZ) 
and consequently every singular vector falls into one of the two rep- 
resentations of that symmetry group. Those labeled I follow the 
identity representation; the omitted columns for Y and Z type ori- 
entations have the same amplitudes as the corresponding (shown) 
columns for X type orientations. The representations labeled 2 are 
twofold degenerate, and take the form A cos <f>x, A cos((f>x + 120°), 
and Acos(cj>x + 240°) for the X, Y, and Z entries, where (A, cj> x ) 
is shown in the table; a second, degenerate singular vector takes the 
same form but with 90° added to all the phases. 



1.5, if the interaction simply scaled as 1/R 3 



2. Singular Value Decomposition 

Applying the singular value decomposition analysis to the 
data in Tables iHl and HI yields the results in Tables HVl and HIT1 
Again, these need to be understood in light of the point sym- 
metries. In particular, every singular vector must transform 
under one of the irreducible group representations, which have 
been identified and indicated in the tables. The irreps of 3 are 
named A (trivial) and E (two-dimensional); the irreps of 2/m 
are named A\ (trivial), A^, B\, and i?2- 

Having identified the largest entries in the interaction ma- 
trix, we can interpret the leading singular values and their sin- 
gular vectors. For example, let's approximate the 6-linkage 
interaction matrix by keeping only its largest terms, found in 
the 4x4 subblock represented by the upper left eight entries in 
TableHB. They are all close to ±17meV, with a pattern of signs 
that gives a rank-1 matrix. There would be just one nonzero 
singular value a = 4(17)meV = 68 meV. Its singular vector 
would have four entries ±1/2, for the four kinds of X orienta- 
tion, and zero on all others. The pattern of signs is (H 1 — ) 

in the convention of Table |IV| which corresponds to the B-2 
representation; in this approximation the top row of the table 
would read (—0.5, 0, 0, 0) and we see this is a good approx- 
imation of the actual singular vector. Similarly, the second 
singular vector in Table II V l is approximately (0, 0, 0.5, —0.5), 
expressing the difference between Z+ and — Z orientations. 

It can also noted that, for the b interaction, both of the two 
leading singular vectors treat the r and I variants orientations 
the same; if we retained only these two terms, that truncated 
interaction sees only only six kinds of orientation, with the r 
and I variants lumped together. On the other hand, the leading 



a (meV) 


irrep 


+X r 


+ Y r 


+ Z r 


-z, 




81.45 




A 2 


0.4621 


-0.1908 










15.59 




A 1 


0.1282 


-0.1625 


0.4881 


-0.4196 




3.836 




B x 


0.2402 


0.3507 


0.2705 


-0.2557 




3.212 




B 2 


0.3727 


0.3333 










2.745 




A 1 


0.3780 


-0.2321 


-0.3247 


0.0329 




1.375 




Si 


0.0484 


0.0077 


0.4245 


0.5612 


2 


733 x 10" 


i 


Ai 


-0.0858 


-0.2940 


0.2701 


0.4894 


2 


073 x 10" 


i 


B 1 


0.3552 


-0.3313 


0.1046 


-0.1312 


1 


364 x 10" 


i 


B 2 


-0.3333 


0.3727 








1 


307 x 10" 


i 


A 2 


0.1908 


0.4622 








2 


017 x 10" 


_> 


Si 


-0.2526 


-0.1308 


0.4854 


-0.3200 




irrep 


signs X 


signs Y 


signs +Z r/l 


signs -Z l/r 








Ai 


+ + ++ 


+ + ++ 


++ 


++ 








B, 


H h 


H h 


+ - 


-+ 








A 2 


+ H 


H h- 


00 


00 








B 2 


H h- 


+ H 


00 


00 



TABLE IV: Singular values a and right singular vectors for b- 
linkage. The singular vectors belong to representations of the point 
group 2m (=C 2v ), as given in the second column: Each singular 
vector has twelve entries, one for each of the orientations (column 
labels). The four orientations beginning with direction "X" or "Y", 
e.g. "(+X r , +Xi, —Xi, —X r )", always have the same amplitudes 
apart from sign, so only the first column is printed from each group 
of four. Similarly, in each of the pairs "(+Z r , —Z r )" and the and 
'+Zi, —Zi)" only one of the two is needed. The pattern of relative 
signs for the group of four or the pair depends on the group represen- 
tation, as given at the bottom of the table. 



singular vector for the c interaction is very far from treating r 
and / equivalently. 



B. Testing validity of discrete effective Hamiltonian 

1. Checks from the fitting process 

We have checked the validity of the orientation potentials in 
multiple ways. First, we report the results of two checks which 
are inherent to the uniform-background scheme for fitting, as 
explained in Sec. HID II 

The first check is that, in the fit with a particular backgroun, 
the apparent one-body energy vector entirely represents inter- 
actions with the background orientation, and therefore has a 
low symmetry. The apparent one-body energy varies between 
+0.03eV and — 0.05eV for most backgrounds. 

Next, from the fitted pair interactions of our chosen clus- 
ters, we can predict the summed pair interactions of either 
cluster with its 13 [100] and [1/2, 1/2, 1/2] neighbors in the 
fixed background, for each of the fixed backgrounds. The 
agreement is excellent: the root-mean-square residual is 2.6 
meV when the varied clusters are a [1 /2, 1/2,1 /2] pair, or 2.3 
eV when they are a [100] pair. This residual is due to either 
multi-body interactions or to more distant neighbor interac- 
tions (with the background clusters). 

Incidentally, in a preliminary study, we directly tested the 
assumption that farther neighbor interactions are significantly 
weaker, by directly fitting V a p trial fits of the pair interaction 
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and therefore could be omitted, for second nearest-neighbors 
displaced by the vectors [0, 1, 1] and also 1, 3]. (This was 
only carried out using a neutral background of Ca-filled clus- 
ters, as mentioned in section HID II ) These interaction ma- 
trices gave singular values more than an order of magnitude 
smaller than those of the nearest-neighbor interactions. 

The second check from fitting is that, insofar as multi- 
cluster interactions can be neglected, different backgrounds 
should give an exactly identical interaction matrix - even if 
there are longer range pair interactions, We found that this 
is largely so. To quantify the dependence on background we 
considered how a particular element V a p of the 12x12 inter- 
action matrix varies while we cycle Q bg through all 12 possi- 
ble values, and compute the standard deviation of V a p. 

In the case of the [1/2,1/2,1/2] cluster linkage, the aver- 
age standard deviation is about 1-2 meV. The larger interac- 
tions and larger standard deviations were between any tetra- 
hedron of orientation +X/Y/Z r and any one of orientation 
-X/Y/Z v 

In the case of the [1,0,0] linkage, the standard deviation was 
typically ~ 0.6 meV), i.e. half as big as for the c-linkage. de- 
viation was larger than average for the 4 2 interaction terms 
that combine X—X type orientations by a factor of almost 
two, and smaller by a factor of nearly two in the 8 2 terms 
not involving any X type orientation. (The X—X pairs had by 
far the strongest interactions, too), 

To summarize, the multi-cluster interactions, as measured 
by the standard deviation of V Q( g while the background is var- 
ied, is ~ 15% of the cluster pair interaction, for both kinds of 
linkage. Thus, the multi-cluster interactions are significant, 
but an order of magnitude down from the dominant pair inter- 
actions. 



the use of constrained relaxation. It should be noted that this 
is entirely canceled out in way we fit the pair interactions (see 
Sec. lnTJ3l 

The scatter of the total energy appears to be a bit under 
±0.1eV, coming from 2 x 3 3 = 54 clusters, or a total of 
54 x 14/2 = 378 linkages, since each has (8 + 6) c- and fe-type 
neighbors. If we imagine that each interaction deviates devi- 
ates randomly and independently from the fitted energy, this 
random energy comes to ~ 4 me V/linkage, consistent in mag- 
nitude with the residuals estimated just above (in Sec. HIIBTI ). 

The fitted Hamiltonian can be further checked for transfer- 
ability to a larger super cell. Figure shows a comparison 
between the actual atomic energies of a random configuration 
in a 4x4x4 with the predictions from the fitted cluster- 
pair Hamiltonian It too shows very high correlation between 
the energies which suggests that the effective Hamiltonian is 
doing a decent job of capturing the essence of the interactions. 



-0.5 



+ + uniform bg 
X X random t-fix 
■ ■ random 
— y=x+c 









-2.0 -1.5 -1.0 -0.5 
EAM Energy(eV)-5.091xl0 3 



2. Tests by direct comparison with total energies 

To fully validate the effective Hamiltonian fitted from uni- 
form backgrounds (in 3x3x3 supercells), we should test 
it on databases with and the energies computed with configu- 
rations other than the fitting database. Our test datasets was a 
3x3x3 supercell with in which randomly oriented tetrahedra 
were placed, and then relaxed in two ways (see Sec. Ill C 2\ , 
giving two variant databsets: 

(1) relaxed with the continuous orientation £2^ held strictly 
fixed; 

(2) the same random orientations, unconstrained relaxation 
so the orientation can change 

(In the latter case, it was assumed without checking that every 
cluster's orientation remains in the discrete well belonging to 
optimal orientation with which it was initialized.) 

Figure[3]shows scatter plots comparing the effective Hamil- 
tonian prediction (vertical axis) with the actual EAM interac- 
tion, for the two random datasets and also the input fitting 
dataset. An offset of about +94.30 eV/cluster has been sub- 
tracted, representing the EAM energy of all the atoms in one 
primitive cell^i A systematic constant offset is visible due to 



FIG. 3: Scatter plot, showing effective Hamiltonian prediction com- 
pared to the relaxed energies, for a 3 x 3 x 3 superlattice with each 
cluster in a random discrete orientation (out of the twelve optimal 
ones). Different symbols refer to different conditions of the fit: + = 
the uniform background configurations (constrained relaxation) used 
for the fit; x = random configurations relaxed with the constraints 
(Sec. Ill C^b and □ = relaxed unconstrained. These last are offset to 
the left, due to the energy reduction when tetrahedra (under forces 
from neighboring tetrahedra) relax to orientations slightly different 
from the discrete minima of the single-tetrahedron orientational po- 
tential. 



C. Cluster Hamiltonian with pair potentials 

All the fits mentioned till now were derived from EAM po- 
tentials. In this section, we have redone them using fitted 
(EOPP) pair potentials. It is not a priori obvious whether pair 
potentials should be adequate for our problem. The Cd2o cage 
atoms make large displacements in response to rotations of the 
tetrahedron, which might have the same cause as the anoma- 
lous c/a ratio in hep elemental cadmium, which (at least is the 
similar element Zn) is poorly captured by the EOPP poten- 
tials^. The response of cage atoms must be the most impor- 
tant determinant of the tetrahedron's elastic interaction with 
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FIG. 4: Scatter plot, comparing energies of relaxed atomic configura- 
tions in 4 x 4 x 4 supercells to those predicted by the discrete cluster 
Hamiltonian, which was fitted from a smaller (3x3x3) supercell. 
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+Y r 
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0.0025 
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0.0044 


0.0025 


— JC r 


-0.0054 
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-Y r 


0.0027 
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0.0017 


0.0044 


— Z r 


-0.0027 


0.0008 
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+x l 


-0.0072 


0.0063 
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+Y l 


-0.0137 


0.0144 


0.0008 


-0.0012 


+Zi 


-0.0158 


0.0144 


-0.0031 


0.0012 


-x t 


0.0126 


-0.0072 


-0.0054 


-0.0055 


-Yi 


0.0220 


-0.0158 


-0.0027 


-0.0021 


-Zi 


0.0220 


-0.0137 


0.0027 


-0.0069 



TABLE V: Interaction matrix for two clusters separated by 
[1/2, 1/2, 1/2] (=c-linkage), using pair potentials (in meV), using 
the same conventions as Table U 



its surroundings. 

The result is that the interactions derived from pair poten- 
tials are remarkably consistent with those from EAM. A sam- 
ple of this is given by the singular value decomposition of the 
c-bond interaction (Table IVDi, in which the top three singular 
vectors show great agreement. 
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-0.0094 
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0.4546 


23.83 
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0.4328 
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0.0286 
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178.0° 
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160.5° 


146.9° 
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E 


0.2354 


0.4099 


0.6458 
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24.2° 
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30.3° 
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-0.2569 


0.4898 
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3.075 x 10~ 2 
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0.2822 


0.4462 


0.2478 


0.5714 
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TABLE VI: Pair-potential results: singular values and right singular 
vectors for c-linkage, in same format as table IIIII 
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0.2241 
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-l 


Si 


0.3844 


-0.2561 
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0.0100 
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-l 


Si 


-0.0314 


-0.2995 


0.4662 


-0.3183 


6.157 x 10" 


-2 


B 2 


0.4259 


0.2619 









TABLE VII: Singular values a and right singular vectors for b- 
linkage, as derived from pair potentials, using the same conventions 
as Table [IV] 



This gives some reassurance as to the independence of our 
results from the specific potentials used. It may help justify 
the adoption of fitted pair potentials for related compositions 
(in particular ScAng) for which EAM potentials are not avail- 
able, subject to the caution that the EOPP potentials must be 
re-fitted if the lattice is compressed or expanded. 

IV. RESULTS: CONTINUOUS ORIENTATIONS 

Till now, the discussion in this article was limited to a set of 
discrete reference orientations, which are determined by the 
local minima of the one-cluster potential. In fact, that is not a 
necessary condition for our analysis: our method of relaxation 
with the rotation constraint (Section lilC 2\ lets us evaluate the 
the effective Hamiltonian for any set of orientations, whether 
stable or not. 

To proceed, we define some discrete set {w a } of m ori- 
entations that samples the continuous space. We adopt the 
uniform background set-up (Sec. HID II letting two interact- 
ing clusters run through all combinations of (w^Wg) while 
keepting the other clusters in fixed orientations. Once again, 
we use the singular-value decomposition (SVD) (12.18b as the 
fitting procedure to obtain the single-cluster and pair terms 
U(£l) and V(Q,£}f). Since there are now thousands of dis- 
tinct pair combinations, the SVD becomes a necessity in the 
case of continuous orientations, rather than an option as it was 
with the discrete version of the interaction. 

The continuous formulation offers multiple opportunities. 
First, we can now map out the single-cluster term, separated 
from the interaction term. This opens up the possibility of de- 
tecting metastable, higher energy minima, of finding the bar- 
riers between discrete wells, and of simulating temperatures 
so high that large deviations from the optimal orientations are 
typical in the ensemble. 

Secondly, we obtain the cluster pair interaction valid for any 
combination of continuously variable orientations, even when 
the interaction term is so strong as to destabilize the local well 
for certain combinations of orientations (or strong enough to 
significantly displace the orientations from the reference ori- 
entation). Finally, going past the SVD this naturally pushes us 
to a further step of simplification in representing the coupling 
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of each cluster, namely orientational harmonics, and these of- 
fer the possibility to unify the basis of functions g(£i) used in 
representing all the different interactions of a cluster. 

In this section, we start off by laying out the quaternion- 
based mathematical framework needed to describe rotations 
(Sec HVBb . Then we carry out two forms of continuous fit. 
First (Sec. lIV Ai l, we limit the rotations of both tetrahedra to a 
single rotation axis. Second (Sec. HVCl l we endeavor to sam- 
ple all rotations. In either analysis, we do find that just the 
first two or three values matter, and the singular vectors 
can be interpolated by smooth functions, thus vindicating the 
motivation of the singular value analysis of the potentials. 



A. "One-dimensional" rotations 

The simplest continuous sub-space of the manifold of 
orientations consists of all rotations around one axis, 
parametrized by a single Euler angle. We choose this axis 
such that the rotations connect (at least) two of the known op- 
timal orientations. Specifically, we choose a zero orientation 
of +X r so that rotation around the x-axis passes through +X r 
to +Xi, —X r , and —Xi. (a rotation by exactly ir takes +X r 
to — X r or +Xi to —Xi.) Rotations around y and z axis are 
also used; we will call the respective rotation angles A x , A y , 
and A z . The sampling points are spaced by 5° along each of 
the three circles. These three data sets are combined into a 
single one, which we call the "three-circle" point set, with a 
total of 36+72+72=180 sample points (the A x rotation runs 
only to 180° since it repeats after that.) 

The first result of the fit (before any SVD analysis) is the 
single-body potential, a byproduct of the processing mandated 
in Sec. HID 31 (Recall that such information cannot be ob- 
tained from the discrete analysis of Sec. [HI]) The results are 
shown in Figure[5la). 

Each deep well is one of the optimal directions. Angle 
A-x.y.z = is, by definition, the orientation +X r . Rotat- 
ing around the x axis, we encounter —Xi at 90°; also, at 
A x w 30° and A — x ~ 120°, we meet +Xi and — X r , respec- 
tively. At either A y = 180° or A z = 180°, we meet —X r ; at 
A y = ±90° we meet ±Z r , while at A z = ±90° we meet 
±Y r . The double well at A x = 105° ± 15° is responsible for 
the ±15° rotation of the tetrahedra in optimal states (of these, 
the well at A x = 90° actually appears to be destabilized by 
the uniform background.) 

In principle, the single-body output has contributions (as in 
the discrete case) from the pair interactions of the background, 
as well as from the true single-cluster term f/(0); however the 
latter contribution is much larger. The background contribu- 
tion merely creates slight energy offsets, visible in the figure, 
between wells which ought to be symmetry-equivalent: e.g. 
along the A^ circle, — X r appears to be lower in energy than 
+X r , by ~ 0.04 eV. (We would eliminate the background 
contribution if we averaged over all possible backgrounds, as 
was done in the discrete analysis of Sec. [nl] but that was not 
carried out in our treatment of the continuous rotations.) 

Next comes the singular vector analysis, with the singu- 
lar vectors normalized according to Eq. (14.3) . The resulting 



^ in frill 

Olll^Ulai ValUC 


c 


- 11 lliva£,C 


U-llllKd^ 


e 




C3 


P600 R30 


C3 P600 


R30 


(Tl 


80.28 


86.04 68.06 


95.20 76.52 


67.86 


0"2 


51.64 


50.87 43.04 


16.17 25.80 


18.40 


03 


28.70 


48.25 31.93 


11.00 6.648 


6.416 



TABLE VIII: First three singular values from three continuous data 
sets (in meV). Here "C3" is the three-circle data set (Sec. IIV At . 
while "P600" and "R30" are the polytope 600 and the random 30 
data sets (Sec. UVCl . 



singular values are included in Table IVIIII The overall mag- 
nitudes are comparable to the discrete result, and some differ- 
ences can be explained away because the three-circle data set 
is not even approximately uniform: e.g., as it starts from +X r , 
it over-represents the four X flavors of discrete orientation. 
(Since those ones have the strongest pair interactions, accord- 
ing to Table Q]] it is not surprising that the dominant singular 
value for the 6-linkage for the C3 data set (Table IVIIIt comes 
out 20% larger than the corresponding singular value for opti- 
mal orientations in Table IIVI) Note also that the third singular 
valus is not so well separated here from the second one, as 
was the case with the discrete orientations (Sec. IIII A 2b . 

In conclusion, the three-circle data set gave a decent indi- 
cation of how smoothly the potential varies between the rel- 
evant discrete orientations, and it goes along with a conve- 
nient way of plotting singular vectors along cuts in this three- 
dimensional parameter space. However, it completely fails 
when we try to transfer to orientations that are not near to 
those three circles in orientation space. Thus, this approach 
does not suffice to provide a parameterized representation of 
the complete functional form for all possible orientations, as 
we would need to use this in a simulation with continuous an- 
gles. 



B. Mathematical handling of continuous rotations 

The most uniform way to parametrize rotations in three di- 
mensions is by four-component unit quaternions: for a rota- 
tion of angle 9 about axis 9, the first component is cos(0/2) 
and the other three are sin(0/2)0. Thus the quaternions map 
out a hypersphere, which corresponds 2-to-l with rotations 
(since antipodal quaternions represent the same rotation). The 
quaternion components are related to the three Euler angles as 
follows: 

go = cos a (4.1a) 
qi = sin a sin 9 cos <fi (4.1b) 
g2 = sin a sin 9 sin 4> (4.1c) 
g3 = sin a cos 9 (4. Id) 

Here 2a is the total rotation angle, and (gi , q2, qs) equals sin a 
times the unit vector of the rotation axis. 

The proper measure in rotation space is uniform on the unit 
hypersphere parametrized by (14. p . Thus, the normalization 
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90 180 270 360 

A, 



FIG. 5: Fitted potentials for a cluster pair separated by [0, 0, 1], dis- 
played along three circular paths in rotation space parameterized by 
rotation angle A, A y , A x respectively, (a), single-body potential term 
U(Ak) along the three circles (the A x plot has period 180°). (b). 
Singular vector <?(£2) for the dominant singular value. (The singular 
vector for the second cluster is related by symmetry to that of the 
first.) Different lines indicate the following databases: (i) Dots [red 
online] = actual measured energies, sampled every 5° for a total of 
180 points; this was the database for the "three-circle" fit; the smooth 
interpolated version of that, using hyperspherical harmonics, closely 
follows the data and is not shown, (ii) solid line [green online] = 
hyperspherical harmonic fits using Polytope 600. (300 orientations 
placed on the vertices of the regular polytope 600 in rotation space, 
fitted by SVD, and then interpolated using hyperspherical harmon- 
ics). A random fit "R150" was practically identical to "P600". (iii) 
The dashed line [blue online] is "R30", a database of 30 randomly 
chosen orientations, symmetrized according to the two mirror planes 
of the [0,0,1] cluster linkage that do not swap clusters, for a total 
database of 120 points. 



convention used for functions of rotation space is 

/t>lTT i>7T i>TT 

d 3 u = { sin 2 ada){ sin OdO) ( / d<j>). (4.2) 
Jo Jo Jo 

In particular, the total volume of cj-space is (2ir) 2 . Hence, we 
normalize a discrete singular vector (gx,g2, •••) by 



m 

^2 ' 



(4.3) 



fc=i' 



If our singular vector is the sampling of a normalized con- 
tinuous mode u(lj) and if the sampling points are uniformly 
distributed, this will be equivalent to the normalization of the 
continuous mode. With such a normalization, singular values 
obtained from different constructions ought to have similar 
magnitudes, as is borne out in Table IVlfll 



1. Hyperspherical harmonics 

We assume the interaction (or singular vector) is a smooth 
function in rotation space. The standard way to parametrize 
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FIG. 6: Same as Figure |5tb), but for clusters related by a "c" type 
separation. Two additional singular vectors are illustrated, as they 
are comparable in strength to the leading ones. Note that this fit is 
less successful, and in places the R30 fit fails, (a) Dominant sin- 
gular vector (b) Second singular vector (c) Third singular vector. 
(The single-body energy plot is very similar to Figure |5}a) and is 
not shown here.) 



such a function with a small number of fitting parameters, is 
a series expansion using some basis of orthogonal functions, 
ideally tailored to the symmetry of the space. When our data, 
sampled at necessarily sparse points in rotation space, is ex- 
pressed in this basis, it can be thought of as simply an elabo- 
rate kind of interpolation. 

The natural basis for the 3-sphere is the hyperspherical har- 
monics, analogous to expanding in spherical harmonics on a 
2-sphere. (These are also commonly used in the theory of 
textures in materials science 2 ^.) We adopt the definitions and 
normalizations for real hyperspherical functions from Eq. (6) 
of Ref. |28l These carry three indices: N, the total "hyperan- 
gular momentum", with L and M to label different functions 
within the same representation. The real hyperspherical har- 
monics, in terms of the three Euler angles a,0,<f>, are then 
given by 



rvMC 

/. 



L (w) 



rvMS 



znlmC^l (cos a)P™ (cos 0) cos M$.4a) 
znlmC^ + _\ (cos a) Pi 1 (cos 0) sin M0t.4b) 



where C^\(cosa) is a Gegenbauer polynomial and 
Pfj (cos 0) is an associated Legendre function. These are or- 
thonormalized with respect to the measure of Eq. (14.2) . and 
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the normalization constant is 



V. ORIENTATIONAL ORDERINGS 



ZNLM = (— 1) 



L+M 



2 L L\ 



(2L + 1) 



(L — M)\ (N + 1){N — L)\ 
{L + M)\ {N + L + 



-,1/2 



(4.5) 



C. Fits from sampling all orientations 

It is necessary to devise some roughly uniform way to sam- 
ple rotation space. We have tried two ways. First, we choose 
m quaternions uniformly spaced by placing them on one of 
the regular polytopes of icosahedral symmetry, the higher- 
dimensional analog of an icosahedron. However, it appears 
the polytope is inefficient because it has too much symme- 
try, and some places in rotation space are far from any point 
of the polytope. Our second approach is to select a random 
list of orientations {u k } and then apply the mm symmetries 
around the [0,0, 1] bond, i.e. those which preserve the two 
cluster positions. This turned out to work much better. 

Applying the SVD yields one dominant singular value, 
along with a corresponding dominant singular vector g^ a as 
shown in Figure [5] (b,c). This represents the leading mode of 
the two-body interaction, and looks like sampling a largely 
sinusoidal function of w Q [Figure [3Jb)]. The y-axis shows a 
noticeably non-sinusoidal profile; 

Fitting the coefficients of hyperspherical harmonic expan- 
sion to the most important singular vector yields, for the 
case of the [0,0,1] (or "b") interaction, coefficients in the 
N = 6 and N = 8 hyperspherical harmonics. Further har- 
monics were not needed and (when included) seemed to re- 
flect sampling arbitrariness and not improve the fit. For the 
[1/2, 1/2, 1/2] (i.e. "c") interaction, we additionally needed 
N = 12 in order to get a decent fit. 

On the other hand, the hyperspherical fits of the 
[1/2,1/2,1/2] ("c") cluster pairs require N = 12 hyper- 
spherical harmonics in addition to the N = 6 and N = 8 
used for the 6-interactions, and still this is not so good a fit. 
Furthermore, the second and third largest singular values are 
non-negligible in the case of the c interaction, as shown in the 
figure. We conjecture that the c interaction is the most com- 
plicated simply because it is the shortest. 

Figure |5jc) shows how the fitted singular vectors compare 
against the singular vector from slices measurements. The 
fits from polytope-600 and the three-circle data set both give 
reasonable approximations of the actual data. In contrast, 
polytope- 1 20 (not shown) fails completely, which can be read- 
ily understood from the fact that the orientations in polytope- 
120 are too spaced too far apart (72° from each other), and 
furthermore the cuts shown in the figures do not even have to 
go through the sampled orientations. 



A great advantage of the fitted Hamiltonian approach is that 
one may discover the optimal structures for systems that were 
not included in the fitting database (and which could not have 
been included, because they are too big). In this spirit, we take 
the discrete effective Hamiltonian from Section [TTT] and find 
what is predicted for the ordering pattern at low temperatures. 
We used Monte Carlo annealing in supercells to discover the 
ground state. 

Watanuki et alA discovered pressure-induced phase transi- 
tions in the CdYbg cubic crystal. Therefore, we extend our 
studies to nonzero pressure, so as to make contact with exper- 
iments that show different ordered states appearing in a range 
of pressures < lOGPa. 



A. Strain, pressure, and bulk modulus 

To extend our calculations to varying hydrostatic pressure, 
we reran calculations of the energies and ground states con- 
strained to various strain values, separated by 0.002 (or by 
0.010 for strains greater than 0.01). The corresponding pres- 
sure P (at T — 0) was then evaluated by recalculating the 
energy at slightly different cell volumes V and using P w 
(AE) I (AV) . The pressure/strain relationship suggests that 
changes in the lattice constant at the T = transitions are 
quite small. 

Note that the "strain 0" results, reported in previous sec- 
tions, used a cell constrained to an a priori lattice constant. 
The present calculation shows that (with the EAM potentials 
we are using) they do not exactly not correspond to zero pres- 
sure, in fact we found zero pressure at negative strain. A pos- 
sible physical meaning to study even more negative strains 
(with negative pressures) is as follows. In isotructural com- 
pounds with the large atom species varied, an increase (de- 
crease) in its effective radius is appears as a negative (positive) 
"chemical pressure", so that we might see a similar phase di- 
agram except for a shift along the P axis. 

The bulk modulus grows rather uniformly with pressure: 
from about 80 GPa to 300 GPa over the range from our most 
negative pressure (strain —0.020) to the largest one (strain 
+0.040). Specifically it was around 150 GPa at P = or 
200 GPa at P « 1 GP. We do not know of any experimental 
measurement of the bulk modulus of CaCdg; measurements of 
the (similar) quasicrystal phase i-CdCa gave a bulk modulus 
68. 1 GPa at zero pressure^. 

We also estimated the pressure using all possible 2x2x1 
supercells with a pair of clusters of all possible orientation 
combinations, placed in uniform backgrounds of all possible 
orientations. This was the same database used to construct 
the cluster pair potentials. These are higher-energy structures, 
so this is roughly like an ensemble of infinite temperature (so 
far as cluster orientations are concerned). We found that (i) 
the pressure was higher by ~ 1.5 GPa, (ii) at zero pressure, 
the strain was —0.008, i.e. an increase of the lattice constant 
which could be interpreted as an orientational contribution to 
the thermal expansion, (iii) the bulk modulus is ~ 20% higher 
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in the ZY structures found at small negative strain, to ~ 10% 
higher around the transition to the ZY structure around strain 
0.07, and roughly unchanged at the highest pressures. 



B. Monte Carlo search for ground state in supercells 

We performed Monte Carlo (MC) simulations using the fit- 
ted discrete cluster Hamiltonian in a 4x4x4 supercell, i.e. 
128 clusters. (Sizes over ~ 50 clusters cannot be equilibrated 
with plain Monte Carlo owing to the frustrated interations.) 
The lowest energy configuration encountered during a run was 
saved and defines the ground state ordering. 



1. Results: ordered states 

The predicted ground state orientation pattern does change 
under pressure. We have found at least three distinct optimal 
states as P is varied. The transitions between them at T = 
are necessarily first-order: since every tetrahedron falls into 
one of twelve discrete orientations, there is no way that one 
pattern can evolve continuously into a distinct one. 

The three T = phases are shown in Figures|7]andr|8] 

(a) Pca2 




(b) C2/c 
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FIG. 7: [Color online] Phases of tetrahedron order in CaCds found 
from simulation at various pressures. [Online: tetrahedra aligned 
with X, Y or Z directions are colored red, green, and blue respec- 
tively; for each direction, a lighter shade is used for both orientations 
in the "+" sense and darker for both in the "— " sense.] They are 
tagged by a provisional label based on the orientations shown in this 
figure, and by their space group (a) "ZY" with 2x2x1 cell, Or- 
thorhombic Pca2 [here in the setting P2ab]. The bcc lattice corner 
clusters have the four directions ±Z r n while body-center clusters 
have the four directions ±Y r n. (b) "Z2" with 2x2x1 cell, Mono- 
clinic C2/c; the 2x2x2 supercell shown is two unit cells. Tetra- 
hedra have two orientations, ±Z r . (c) "XYZ" with 4x2x2 cell, 
CHECK space group in the setting ... The structure is shown as two 
slabs, each one lattice constant thick. 

We saw three patterns in the whole range of pressures (ac- 
tually strains): 

(1) at strain -0.020, i.e. roughly P < -2.3GPa, we saw 
a structure we call "ZY" with a 2x2x1 unit, thus 
8 clusters per cell, which appears to have orthorhombic 
space group Peal [Figure |7ja)] 



(2) for strains -0.018 to -0.004, i.e. pressures -2.3 to 
—0.6 GPa, we found a simpler structure we call Z2 con- 
taining just two cluster orientations which simply form 
(110) layers. [Figure |7Jb)]. The space group is cen- 
tered monoclinic C2/c with 4 clusters per cell (=2 per 
primitive cell). 

(3) For strains —0.020 through 0.008, i.e. up to a pres- 
sure of ~ 1.4GPa, we see the complex 4x4x2 
pattern shown in [Figure |7jc)], which has a triclinic 
(pseudo monoclinic) unit cell; there is a centering oper- 
ation within the 4x4x2 cell so there are 32 clusters 
per cell. 

Finally, for the large strain of 0.010 (around P = 1.6 GPa), 
we found a less regular 4 x 4 x 2 ; we do not know if this is 
the true ground state. 

We will describe the 4x3x2 ordering pattern as a stacking 
of 4 x 4 layers along the z direction. The structure repeats 
after four of these layers (i.e. after two lattice constants, as 
there are layers of cell corners alternating with layers of body- 
center sites). We will use "flavor" to designate the X, Y, or 
Z nature of the orientation, so that each flavor includes four 
orientations. The even layers (cell coiners) are of one kind 
that we call "XY" layers, after the orientation flavors found 
in them. The odd layers (body centers) are of another kind we 
will call "XZ" layers. 

The basic ingredient of a layer is a chain with a pattern of 
orientations AAAA where A stands for inversion of A, where 
A = X, Y, or Z. Every layer is made by stacking such four 
chains side-by-side, using two alternating flavors; the second 
occurrence of the same flavor uses the orientations not found 
in the first one. In the XZ layers, the chains run in the x direc- 
tion and are stacked in the y direction. Note that, if we draw 
bonds between neighbors of identical orientations, this forms 
a columnar pattern of dimers covering the square lattice. The 
XY layers have chains the other way around, running in the y 
direction and stacked in the x direction, but slightly different 
from the way in the XZ layers. The slight difference is that 
here, if we mark the adjacent pairs of the same orientation, it 
forms a staggered dimer pattern. 

If we go up two layers, i.e. one lattice constant in the z 
direction, we get the same pattern, except all orientations are 
replaced by the complementary ones of the same flavor. Also, 
the registry between successive layers is such that (in projec- 
tion) the X dimers in the XZ layer cross the Y dimers in the 
XY layer. 

We also did simulations with smaller dimensions of super- 
cell. Runs with a 2x2x1 cell indeed give the same 
ground state as the 4x4x4 whenever that state fits into 
the smaller supercell, i.e. for negative strains; for positive 
strains, the 2x2x1 gives higher energy states, showing 
that the ground state cannot fit into that cell. (Monte Carlo 
results from a 3x3x3 system are always higher in energy, 
presumably because the true ground state can never fit into an 
odd dimensioned cell.) Finally, we also exhaustively enumer- 
ated all 12 s orientation combinations in the 2x2x1 super 
cell. The ground states found in this way agree with the MC 
results whenever that shows a ground state with a 2x2x1 
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FIG. 8: The same structures shown in Figure [7] expressed in terms 
of the orientation labels defined in Sec, III ATI 



supercell, i.e. at strains less than or equal to zero. 

Brommer performed an exhaustive enumeration using a 
x 2 su perc ell, (this contains 8 clusters). His best 
structure (see Ref.[l9l Figure 4.13) appears to be our a phase, 
containing the four orientations ±Z r n, which is similar but 
not identical to our result at zero pressure, which was the "Z2" 
phase. That is surprising, in view of the close similarity of our 
fitted potential to that of Brommer et al (see SecfA] above). 



2. Comparison with experiments 

Our results are reminiscent of, but not in agreement with, 
the experiments, which find ordered states with unit cells 
of either 2x2x1 or 2x2x2. The experimental 
pressure-temperature (P-T) phase diagram of Watanuki et 
alii is shown in Figure [9] The ground state structure shown 
in Figure |7tb) agrees with the C2/c proposed arrangement in 
Ref.|22| similar to one in Ref. [3 

Direct comparison to experimental data is not possible 
since the P > experiment is for a different alloy CdgYb. 
Still, some similarities exist. Most importantly, the Z2 phase, 
which we predict to be stable in the pressure range bracket- 
ing P = 0, is (apart from small rotations around z) the same 
C2 / c structure proposed in Figure 5 of Ref. |30|for the (ambi- 
ent pressure) phase of CaCd6- The {110} structure suggested 
in Figure 3(c) of Ref. 3, for the low- and high-pressure phases 
of CdgYb, is also the same as our or Z^ phase if we use 
zero rotations in place of the small right/left rotations. 

We have not studied our model system at the highest pres- 
sures, so we do not know whether it has a counterpart for the 
highest-pressure phase of Ref. 0. 

Watanuki et at suggest that the reason that different phases 
are selected with increasing pressure is the competition of 
nearest-neighbor and longer-range interatomic interactions. 
The latter were due, they conjectured, to Friedel oscillations. 
(We point out that elastically mediated interactions are just 
as long -ranged as the oscillating part of the interatomic po- 
tential, both of these being ideally l/R 3 .) However, we ob- 
tained a similarly complex phase diagram without refitting the 
EAM potentials (meaning that any change in Friedel oscilla- 
tions was not explicitly taken into account). Within our ap- 
proach, the reason for the multiple phases is that the orien- 
tational interaction of neighboring clusters is quite frustrated: 
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FIG. 9: (a) Pressure-temperature phase diagram of CdsYb [after 
Watanuki et aP, Figure 3(a)]. The ordering wavevector of each phase 
is indicated (higher temperature phases have similar order but par- 
tially disordered); region with {111} phases is shaded. All transi- 
tions to the disordered (bcc) phase are continuous, (b) Our predic- 
tion for CaCds; note shift in the pressure axes. Critical pressures at 
T — are indicated by vertical bars; critical temperatures with first- 
order transitions, as reported in Table UXI for the three strains where 
this was measured, are shown by black circles. 



low-energy ordering patterns cannot simultaneously satisfy all 
interactions. It seems there are several inequivalent ways to 
balance good interactions with bad interactions that are nearly 
degenerate. Thus, a small change in the relative cost of differ- 
ent orientation combinations is expected to tip the balance to 
a different pattern^ 



C. Transitions at T > 

Tamura et found an order-disorder transition in 

CaCdgat T m 100K. Under pressure, in the similar system 
CdgYb, Watanuki et al 4 found order-disorder transitions of 
the high-pressure phases too; order set in at T =200 - 250 K, 
and a further ordering within that phase occurred at T « 140 
K. 

Simulations by Brommer et al based on their version of the 
discrete cluster pair Hamiltonian, using a 4x4x4 simu- 
lation cell, found a first-order transition at T w 9 IK. Brom- 
mer pointed ouli2 this represents an entropy jump of about 
1 .0 ks per cluster - twice what was estimated by Tamura et 
al^ for the similar compound Cd@ Y - and an energy jumo of 
~ lOmeV/cluster. 

At present we only have preliminary data concerning tran- 
sitions as a function of temperature. The tempering MC sim- 
ulation, which requires running multiple replicas of the sys- 
tem at different temperatures, naturally detects discontinuities 
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strain 


P 


T c 


AE 


AS 




(GPa) 


(K) 


(meV/cluster) (k 


B/cluster) 


-0.010 


-1.46 


72.5 


3.359 


0.537 






97.5 


6.953 


0.827 


0.000 


-0.14 


40.0 


0.547 


0.159 






70.0 


3.672 


0.609 


+0.010 


1.65 


100.0 


6.641 


0.770 



TABLE IX: Thermal transitions at three fixed strains, along with the 
corresponding pressure at T = 0. The T c is uncertain by ±2.5 K. 
The jump in total energy is shown for each transition, divided by the 
number of clusters, and the corresponding entropy change AS per 
cluster is inferred (in units of Boltzmann's constant). 



in the energy as a function of temperature; our present re- 
sults, shown in Table llXI . are based purely on this metric. We 
took data from zero to high temperatures for three different 
choices of the strain: strongly negative, zero, and strongly 
positive, finding two first order transitions for the first two 
cases, but only one in the case of positive strains. Around 
zero strain, there seems to be a particular tendency to have 
closely competing states and low-energy excitations, and we 
believe this explains the one rather low transition temperature 
for that case. 

Simulations by Brommer et al^- furthermore found a ther- 
mal ordering transition in a 4x4x4 supercell; they did 
not analyze the orientation pattern in the ordered state there, 
but found its energy at the transition was only 1 meV/cluster 
higher than the s/2 x s/2 x 1 structure. 

Presumably, both in experiment and in our simulations, the 
high-T phases are partially disordered versions of the low- 
T phases seen at the same pressure. In particular, the ori- 
entations related by a change of r to I suffix in their sym- 
bols differ by a comparatively small rotation, so possibly (see 
Sec. lIII A~TT i they have similar interactions. Thus one possibil- 
ity is that a partially disordered states has a sublattice contain- 
ing (say) random, equal populations of +X r and +Xi orien- 
tations. 

In many ternary cases, also in ScZn 6 , the ordering tran- 
sition is not sharp in experiments, which Tamura speculated 
indicates a glassy freezing. This is plausible in view of the 
frustrated orientational interactions we found. 



VI. DISCUSSION 

We will first summarize this work and then the outlook for 
extensions of it. 



1. Summary 

In this paper, we presented a comprehensive template for 
studying cluster interactions in CaCdg and more generally in 
any material possessing cluster orientation degrees of free- 
dom. We showed how to operationally relate cluster orien- 
tations to atom positions (Sec. Ill Cb . worked out some group 
theory for the symmetry of the interactions (Sections HID 21 



and IIII Al l, and fixed the technical obstacle of redundant pa- 
rameters due to constraints in the pair counts (Sec. HID 3) . 
Most of all, we showed that a singular value decomposition 
can clarify the dominant nature of the interaction and allows a 
long list of fitting parameters to be truncated to a manageable 
number (Sections III El andl HI A 2\ . 

When this was actually applied to CaCdg, we found (in 
Sec. IIII B I that the omitted interactions - of whatever form, 
multi-cluster or pair interactions with farther neighbors - 
amounted to only 1/100-1/10 of the nearest-neighbor pair in- 
teractions. The same is true even for pair interaction contribu- 
tions apart from the first two singular vectors. This means that 
truncating to those two dominant terms, as about as realistic 
as the nearest-neighbor Heisenberg spin exchange, typically 
is, when used to model a magnetic material. 

With effective interactions in hand, we carried out Monte 
Carlo simulations of larger lattices, for the purpose of discov- 
ering the optimum arrangements (Section|V}. This is difficult, 
as the interactions are frustrated leading to glassiness. Along 
with this, we made a rather sketchy study of the phase dia- 
gram, varying pressure and temperature. Its overall nature is 
broadly reminiscent of experiment (on CaCdg or isostructural 
compounds such as CdgYb), in that several ordering patterns 
are encountered as pressure is varied, and there are also multi- 
ple transitions with increasing temperature, which we respect 
represent partial orderings. 

Using our method of orientation-constrained relaxation 
(Sec. Ill C 21 ), it was possible to extend our analysis to contin- 
uously varying orientations of the clusters. For that case, aug- 
ment couple the singular-value-decomposition with a decom- 
position in terms of rotational harmonics. We have not fully 
developed all that could be done with continuous orientations. 
One byproduct of this calculation, which was not available 
otherwise, is the single-cluster orientational potential energy. 
If the continuous type interactions can be parametrized use- 
fully for Monte Carlo simulations, one application would be 
to study the dynamics of clusters including the paths by which 
they pass from one discrete well to another. (Such a study was 
done for ScZng using all-atom molecular dynamics with pair 
potentials^.) 



2. Possible future work 

One obvious direction for future work is to study the com- 
position dependence. For example, diffraction found some- 
what different ordering patterns as one varies the large atom 
component in isostructural alloys (Ca, Yb, and the other rare 
earths have slightly different sizes): e.g. s/2 x s/2 x 2 cell in 
1/1-CdgYb, a versus s/2 x s/2 x 1 C-centered monoclinic cell 
in l/l-CdgGa. Similarly, different critical temperatures were 
measured experimentally. Certainly, the strength of tetrahe- 
dron interactions will depend sensitively on the composition. 
In ScZn(, 12 the interactions are much weaker than in CaCdg 
but an ordering occurs nevertheless 10 . 

Another direction is to thoroughly study the thermal behav- 
ior, which will require metrics to identify the nature of par- 
tially ordered states. This can also be extended to "approxi- 
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mant" phases with larger cells such as the 2/1 approximant of 
i-CdYb, which has four equivalent clusters per cubic cell; this 
shows a transition to tetrahedron orientational order, a com- 
plex stacking along a (100) directional. 

Our ultimate goal is to understand the tetrahedron orienta- 
tions in the quasicrystal phase, and their role in stabilizing it. 
A speculative possibility is the implementation of matching 
rules by such clusters. Matching rules in the Penrose tiling 
(or its 3D analog) are markings that spoil the symmetry of 
otherwise rhombic or pentagonal objects, and enforce a de- 
terministic, quasi-periodic arrangement analogous to an ideal 
crystal, thus offering one scenario for the stabilization of qua- 
sicrystals. It should be noted that the atomic structures of 
various icosahedral quasicrystals, understood as packings of 
fully symmetric clusters, showed no features that could im- 
plement such matching rules. Furthermore, there is a more 
economical alternative scenario that long-range order is emer- 
gent in dimension 3 from a "random tiling" in which the clus- 
ter packings sample an ensemble of extensive entropy. The 
random-tiling scenario found support (in various icosahedral 
quasicrystals, including i-CaCd) in the shapes of diffuse tails 
around Bragg peaks in diffraction experiments . Otherwise 
no decisive experiments are known, so the best approach to 
discover matching rules, if they exist, appears to be multi- 
scale simulation of the sort carried out in the present paper. 
The most plausible specific physical mechanism for matching 
rules in an icosahedral quasicrystal is via asymmetric inner 
clusters such as the tetrahedron in i-CaCd. These might be in- 
vestigated in the future by extending the methods of this study 
to larger approximants such as Cai3Cd 7 g2L. 

Acknowledgments 

This work was supported by DOE Grant DE-FG02-89ER- 
45405 (WC, CLH, MM), and Slovak Research and Develop- 
ment Agency funding under contracts VEGA 2/0111/11 and 
APVV-0076-1 1 (MM). We are grateful to P. Brommer for pro- 
viding the interaction matrices computed in Refs.[l6landll9l 

Appendix A: Comparison to Brommer et al fit 

We were provided the interaction matrices computed by 
Brommer et al&. They had recognized the redundancies in 
the fit, and resolved them arbitrarily by setting certain terms 
to zero. That meant the magnitude of any particular pair inter- 
action V a p could not be given a physical interpretation: such 
interactions are valid only for computing the total energy of 
an entire configuration. However, if one applies our recipe to 
resolve redundancies from Sec. HID 3l to their V a p, we get a 
well-defined result which can be compared with ours. 

We expect a great similarity to our EAM results, since we 



used the same potentials as they did. The important differ- 
ences in our calculation are 

(a) For a database, out of the approaches we described 
in Sec. HID 11 they used random whole configurations 



whereas we used the uniform background. 



a 


irrep 


+X r 


-Xi 


+x t 


— JC r 


27.69 


E 


0.5542 


0.5542 


0.1617 


0.1617 






3.9 


-3.9 


-160.6 


160.6 


32.62 


E 


0.3927 


0.3927 


0.4232 


0.4232 






-91.2 


91.2 


125.8 


-125.8 


22.6 


E 


0.4139 


0.4139 


0.4025 


0.4025 






-104.3 


104.3 


-148.9 


148.9 


18.25 


I 


0.3401 


-0.3401 


-0.2259 


0.2259 


13.15 


I 


0.2887 


0.2887 


-0.2887 


-0.2887 


7.385 


E 


0.1842 


0.1842 


0.5472 


0.5472 






44.8 


-44.8 


48.8 


-48.8 


5.404 


I 


0.2259 


-0.2259 


0.3401 


-0.3401 



TABLE X: Results from Brommer et al, Ref. [If]: singular values 
and right singular vectors for c-linkage, same format as table [III] 
Note the symmetries relating the pair of columns for +X r and — X r , 
and similarly relating the columns for +Xi and —Xi : they have the 
same magnitudes and, for representation E, they have opposite phase 
angles (after removal of an overall arbitrary phase). 



(b) They did not use constrained relaxation, but relied on 
the final states corresponding to the initial ones. (See 
our discussion in Sec. HIC2l ) Furthermore, their initial 
state was not one of the twelve idealized orientations 
inferred from relaxations, but one of the Gomez-Lidin 
orientations, postulated on the basis of diffraction. 

The results look similar for the 6-linkage interaction ma- 
trix. When we perform the SVD analysis, the overall pattern 
of singular values seen in Table IXll is similar to what we saw 
in Table llVb . though with slightly less of a separation between 
large and small singular values. However, the dominant sin- 
gular vector comes from a different representation, and the 
pattern of signs in the singular vectors is also different. (It is 
likely these differences are artifacts of different conventions 
we used for the orientations; Refji& used the other choice for 
orienting icosahedral axes in cubic symmetry, differing from 
ours by a 90° rotation. 

On the other hand, the results for the c-linkage interaction 
matrix look completely different, as does the singular value 
decomposition (Table [X]>. Furthermore, for reasons we do not 
understand, the latter exhibits an additional symmetry relat- 
ing the columns differing by the sense + and — of the ori- 
entation, which does not correspond to any symmetry of the 
clusters. Finally, the c-linkage singular values are large and 
(unlike all other cases analyzed in this paper) the magnitudes 
of the largest and smallest ones do not differ by a large factor. 
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a (meV) 


irrep 




+Y r 


+ Z r 


-Zi 


69.37 


B 2 


-0.1587 


-0.4741 








21.75 


Ai 


-0.0851 


0.1515 


0.3978 


-0.5305 


6.194 


Ax 


0.0064 


0.3327 


-0.4950 


-0.1832 


5.415 


B 2 


-0.4741 


0.1587 








4.787 


Si 


0.4921 


0.0885 








4.280 


A 2 


0.1681 


0.1638 


0.4602 


-0.4220 


2.415 


A 2 


0.0087 


-0.4041 


0.3967 


0.1259 


1.672 


A 2 


0.1684 


-0.2417 


-0.3588 


-0.4446 


8.392X1CT 1 


Si 


0.0885 


-0.4921 








7.456xl0 _1 


A 2 


-0.4396 


-0.0380 


0.0464 


-0.3292 


7.442xlCT 2 


Ax 


0.3992 


-0.1818 


-0.1160 


-0.3189 



TABLE XI: Results from Brommer et al, Ref. [3: singular values 
and right singular vectors for fa-linkage, same format as table llVI 



A.P. Tsai, J. Q. Guo, E. Abe, H. Takakura, and T. J. Sato, Nature 
408, 537-8 (2000). 

H. Takakura, J. Q. Guo, and A.-P. Tsai, Phil. Mag. Lett. 81, 411 
(2001). 

J. Q. Guo, E. Abe, and A. P. Tsai, Phys Rev. B 62, R14605 (2000). 
T. Watanuki, A.Machida, T.Ikeda, K.Aoki, H.Kaneko, T.Shobu, 
T. J. Sato, and A. P. Tsai, Phys. Rev. Lett. 96, 105702 (2006). 
R. Tamura, "Properties of Cd-Based Binary Quasicrystals and 
Their 1/1 Approximants" Isr. J. Chem. 51 SI 1263-1274 (20 11). 
M. Mihalkovic and C. L. Henley, preprint larXiv: 1 1 12.38041 . 
K. Nozawa and Y. Ishii, J. Phys. Condens. Mat. 20, 315206 
(2008). 

"Experimental evidence for a phase transition inaZneSc 1/1 cubic 
approximant" T. Yamada, R. Tamura, Y. Muro, et al. Phys. Rev. B 
82, 134121 (2010). 

K. Nishimoto, R. Tamura, and S. Takeuchi, Phys. Rev. B 81, 
184201 (2010). Phys. Rev. B 81, 184201 (2010) [5 pp] "In situ 
transmission electron microscopy observation of an orientational 
order-disorder transition in CdeEu and CdsCe crystalline approx- 
imants", 

H. Euchner, T. Yamada, H. Schober, S. Rols, M. Mihalkovic, R. 
Tamura, T. Ishimasa, and M. de Boissieu, J. Phys. Condens. Mat- 
ter 24, 415403 (2012), "Ordering and dynamics of the central 
tetrahedron in the 1/1 ZnsSc periodic approximant to quasicrys- 
tal" 

T. Hatakeyama, K. Nozawa and Y. Ishii, Z. Kristallogr. 223, 830- 
832 (2008) 

M.Mihalkovic and C. L. Henley, Philos. Mag. xx, 1 (2011). 
C. P. Gomez and S.Lidin, Phys. Rev. B. 68, 024203 (2003). 
A. Palenzona, J. Less-Common Met. 25, 367 (1971). 
P. Brommer and F. Gahler, Philos. Mag. 86, 753-758 (2006). 
P.Brommer, F.Gahler, and M.Mihalkovic, Philos. Mag. 87, 2671 
(2007). P. Brommer, F. Gahler, and M. Mihalkovic, Philos. Mag. 
87, 2671-2677 (2007). 

Marek Mihalkovic and C. L. Henley, Phys. Rev. B 85, 092102 
(2012). 

S-Y Piao, C. P. Gomez, and S. Lidin, Z. Naturforsch. 60b, 644-9 
(2006,) 



P. Brommer, Ph. D. thesis (Universitat Stuttgart, 2008): "Devel- 
opment and Test of Interaction Potentials for Complex Metallic 
Alloys" 

E. R. Grannan, M. Randeria, and J. P. Sethna, Phys. Rev. Lett. 60, 
1402 (1988). 

C. P. Gomez, Angew. Chemie 40, 4037 (2001). 
R. Tamura, K. Nishimoto, S. Takeuchi, K. Edagawa, M. Isobe, 
and Y. Ueda, Phys. Rev. B 71, 092203 (2005). 
J. R. Shewchuk, "An Introduction to the Conjugate Gradi- 
ent Method Without the Agonizing Pain", online resource 
http://www.cs.cmu.edu/quake-papers/painless-conjugate-gradient.pdf 
page 52. 

C. L. Henley, Phys. Rev. B43, 993 (1991). 

J. Richmond-Decker, B.S. thesis (Cornell University, 2012); J. 
Richmond-Decker, D. Chaudhuri, M. Mihalkovic, and C. L. Hen- 
ley, unpublished. 

Uzi Hizi, Ph. D. thesis, Cornell University (2006). 

M. Mihalkovic, J. Richmond-Decker, C. L. Henley, and M. 

Oxborrow, "Ab-initio tiling and atomic structure for decagonal 

ZnMgDy quasicrystal" (preprint 2011) |arXiv:1112.395l| 

J. K. Mason and C. A. Schuh, Acta Materialia 56, 6141 (2008). 

J. Z. Jiang, L. Gerward, and J. S. Olsen, Appl. Phys. Lett. 79, 2538 

(2001). 

R. Tamura, K. Edagawa, K. Shibata, K. Nishimoto, S. Takeuchi, 
K. Saitoh, M. Isobe, and Y. Ueda, Phys. Rev. B 72, 17421 1 (2005). 
R. Tamura, K. Minoda, S. Takeuchi, K. Edagawa, T. Kiss, T. 
Yokoya, and S. Shin, Philos. Mag. 86, 489-497 (2006). 
R. Tamura, Y. Murao, S. Takeuchi, M. Ichihara, M. Isobe and Y. 
Ueda, Jap. J. Appl. Phys. 40, L912 (2002). 

R. Tamura, K. Edagawa, C. Aoki, S. Takeuchi, and K. Suzuki,: 
Phys. Rev. B, 68, 174105 (2003). 

Our definition is intended strictly for T = 0. At finite temperature, 
the effective Hamlitonian should properly be defined as a free en- 
ergy, including both energy and entropy contributions, and sam- 
pling all configurations with the given tetrahedron orientations. 
That would obviously be much harder to evaluate from simula- 
tions. 

The "classical" potentials are of course by anchored fitting to 
electronic calculations that are based, via the Local Density Ap- 
proximation, to the many-electron Schrodinger equation. They are 
classical in the sense that they are deterministic functions of the 
atomic positions, having integrated out the electron wavefunctions 
and neglected the ionic zero-point motions. 
Actually, our resolution of the redundancy is equivalent to pro- 
jecting away the subspace in V a p parallel to (1, 1, 1, 1...) so its 
rank was reduced by one. Therefore, the last singular value a m is 
always zero. 

Interestingly, Figure[3]demonstrates that any uniform background 
state is about 20 meV/cluster higher in energy than a typical ran- 
dom pattern. 

A change in the relative strength of (1/2,1/2,1/2) versus 
(1, 0, 0) type cluster pair interactions could be part of this change. 
This might be similar in its mechanism to the change posited in 
Ref. 0, but is not the same, since they imagined a change in the 
ratio of hard-core and long-range interactions. 



